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Abstract 

The full non-linear evolution of the tidal instability is studied numerically in an ellipsoidal fluid domain relevant for 
planetary cores applications. Our numerical model, based on a finite element method, is first validated by reproducing 
some known analytical results. This model is then used to address open questions that were up to now inaccessible 
using theoretical and experimental approaches. Growth rates and mode selection of the instability are systematically 
studied as a function of the aspect ratio of the ellipsoid and as a function of the inclination of the rotation axis 
compared to the deformation plane. We also quantify the saturation amplitude of the flow driven by the instability and 
calculate the viscous dissipation that it causes. This tidal dissipation can be of major importance for some geophysical 
situations and we thus derive general scaling laws which are applied to typical planetary cores. 
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1. Introduction 



Tides are large scale deformations induced by gravitational interactions that affect all the layers of a given planet or 
star. They play an important role in geo- and astrophysics and so, many studies are devoted to this subject. Their most 
obvious phenomena are of course the oceanic flows on Earth, but tides are also responsible for the intense volcanism on 
Io for instance. In stars but also in liquid planetary cores, tidal forcing induces an elliptical deformation of the rotating 
streamlines that may excite a parametric resonance of inertiai waves called the elliptical or tidal instability. This 
instability could, for instance, be responsible for the surprising magnetic field in Io llzj, l26l [l3ll and for fluctuations 
in the Earth's magnetic field on a typical timescale of 10,000 years IjJ]. It may also have a significant influence on the 
orbital evolution of binary stars Irnll and moon-planet systems l3~lll . 

As described in the review of |22|], the elliptical instability is a three-dimensional instability which may grow 
as soon as a rotating fluid possesses elliptical streamlines. As previously stated, it is due to a triadic parametric 
resonance between two inertiai waves of the rotating fluid 112411 and the underlying strain field responsible for the 
elliptic deformation ylH]]]. Such an instability has been found in many different contexts, including those where 
the strain field is due to vortex interactions or to elliptically deformed boundaries. It has thus been the focus of 
numerous theoretical and experimental studies, devoted to the dynamics of two-dimensional turbulent flows (e.£ 
162, 45 1 or lfl7ll for a viscoelastic fluid), to the stabili ty of wakes 
and to the flow inside deformed rotating cylinders 
recently, magnetohydrodynamical effects of the elliptical instability have been investigated, with applications in MHD 
turbulence [58 



icvoieu lo me uynamics oi iwo-umieiisionai luiuuieiu nows l^e.g. 
ibility of wakes ll37ll . to th e dy namics of vortex pair s ll38i[35l43n . 
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or in induction processes in planets II26L 1 1 311 . 



From a numerical point of view, most studies have been devoted to the dynamics of deformed two-dimensional 
vortices in three-dimensional domains [e.g. (39[ [56l |29l 52]. For instance, the seminal paper of [48] considered the 
growth of a perturbation on the planar velocity field associated with an elliptical vortex inside a box with zero normal- 
flow boundary conditions, assuming periodicity in the axial direction; the linear eigenvalue problem was then solved 
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by spectral methods, and the growth rate of the elliptical instability was found. Numerical studies in rotating containers 
are less numerous. 14211 used a non-orthogonal elliptico-polar coordinate system in order to solve the flow in an 
elliptically deformed cylinder by spectral methods; the non-linear temporal evolution of two different modes of the 
elliptical instability was calculated with no-slip boundary conditions on the sidewalls and stress-free conditions on 
the top and the bottom of the cylinder. Il54ll determined numerically the frequencies and the growth rates of the 
instability in both ellipsoid and ellipsoidal shell geometries using a linear Galerkin method, projecting the flow on a 
selected number of inertial waves. Finally, [55] studied the stability of self-gravitating compressible ellipsoidal fluid 
configurations and pointed out the occurence of the elliptical instability. However, to the best of our knowledge, the 
full non-linear evolution of the tidal instability in an ellipsoidal geometry has not yet been simulated numerically. This 
is th e purpo se of the present work, which aims at completing our previous theoretical and experimental investigations 
0271 l26l 1 13l 13111 in accessing global quantities of fundamental importance for planetary applications. 

In this paper, we focus on hydrodynamics only, leaving thermal and magnetic field interactions with the tidal 
instability for further investigations [e.g.@]. In section [2] the system and the numerical method are presented. The 
numerical model is validated by comparison with the theoretical results found in the literature, regarding the flow in a 
rotating sphere submitted to a small and fixed elliptical deformation. In section[3] we then study the influence of three 
additional complications of fundamental importance for geophysical applications: (i) the influence of the length of the 
ellipsoid along its rotation axis compared to the mean equatorial radius (i.e. the body oblateness in geophysical terms), 
(ii) the influence of a background rotation, corresponding to the orbital motion of the companion body responsible 
for the tidal deformation, and (iii) the influence of the inclination of the rotation axis compared to the deformation 
plane (i.e. the obliquity). In section [4] we quantify systematically two large-scale quantities relevant for planetary 
applications, namely (i) the amplitude of the instability at saturation, and (ii) the power dissipated by the instability. 
General laws in terms of dimensionless numbers are derived and applied to typical planetary cores in section|5] 



2. Numerical model and its validation 

2.1. Definition of the system and description of the numerical method 

The present study takes place in direct continuity of our experimental studies of the elliptical instability in a 



deformed spheroid 1271 1261I13LI3 111 . In these experiments, a hollow sphere of radius R, molded in a silicone cylinder, is 
filled with liquid and set in rotation at a constant angular velocity D. about its axis (Oz), while it is slightly compressed 
by a quantity s along the axis (Ox), perpendicular to the rotation axis. The geometry is then a triaxial ellipsoid of axes 
(a, b, c) — (R— s, R + s, R), with an equatorial ellipticity e = (b 2 - a 2 ) / (a 2 + b 2 ) and a constant tangential velocity along 
the deformed boundaries, equal to OR at the equator. Such a configuration is a model for a liquid planetary core with 
no solid inner core (e.g. like the jovian moon Io or the early Earth), surrounded by a deformable mantle with a constant 
tangential velocity. Similarly, our numerical model studies the rotating flow inside an ellipsoid of axes [a, b, c) related 
to the frame [Ox, Oy, Oz), with a constant tangential velocity all along the boundary in each plane perpendicular to 
the rotation vector £2 (see the sketch in figure[T]l. In order to extend our previous experimental study, the length c of 
the ellipsoid can be chosen independently of the other lengths a and b (with b > a). Moreover, the rotation axis of 
the ellipsoid can be inclined compared to the c-axis. Note however that in this paper, unless otherwise specified, we 
choose c equal to the mean equatorial radius R eq — {a + b)/2 and a rotation axis along (Oz), as in the experimental 
setup. In all simulations, the fluid is initially at rest and we suddenly impose at time t — a constant angular rate 
Q. such that the tangential velocity along the deformed boundaries in each plane perpendicular to the rotation axis, is 
equal to Q. (a' + b')/2, where a' and b' are the axes of the elliptic boundary in this plane. In the following, results 
are non-dimensionalised using the mean equatorial radius R eq as the length scale and Q _I as the time scale. Then 
five dimensionless numbers are used to fully describe the system: the Ekman number E — v/(Q R 2 eq ), where v is the 
kinematic viscosity of the fluid, the ellipticity e — (b 2 - a 2 ) /(a 2 + b 2 ) of the elliptical deformation, the aspect ratio 
c/b which quantifies the oblateness of the ellipsoid, and finally the inclination 9 and declination cp of the rotation axis. 
The problem numerically solved is then described by the following system of dimensionless equations: 

— + u Vu = -Vp + E A u — 2 £2* x u, (1) 

at 



V • u = 0. 
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(2) 



where the no-slip boundary conditions are used for the fluid. Note that we work is the reference frame where the 
ellipsoidal shape is at rest, which is the inertial frame of reference in most of the paper. The Coriolis force -2 fi'XH 
is only used in section [3721 where the whole triaxial ellipsoid is submitted to a global rotation at Q.* e z . 

Usually, numerical studies of planetary cores benefit from their spherical geometry to use fast and precise spectral 
methods. In our case however, there is no simple symmetry. Our computations are thus performed with a standard 
finite element method widely used in engineering studies, which allows to deal with complex geometries, such as 
our triaxial ellipsoid, and to simply impose the boundary conditions. Note that a very efficient finite element method 
was recently introduced by J9[], but it is up to now restricted to spheroidal geometries (a = b). Using a commercial 
softwar^H, an unstructured mesh with tetrahedral elements was created. The mesh element type employed is the 
standard Lagrange element PI - P2, which is linear for the pressure field but quadratic for the velocity field. Note that 
no stabilization techniques have been used in this work. We use the so-called Implicit Differential-Algebraic solver 
(IDA solver), based on backward differencing formulas 111 811 . At each time step the system is solved with the sparse 
direct linear solver PARDISC@. 

The elliptical instability induces a three-dimensional destabilization of the initial two-dimensional elliptical stream- 
lines. To study its global properties, it thus seems natural to introduce the mean value of the vertical velocity 
W = ^ JJJ V \w\ dr, with w the dimensionless vertical velocity and V the volume of the ellipsoid. The typical evolution 
of W as a function of time is shown in figure|2](a) for E = 1/500 and s = 0.317. At t = 0, the fluid is at rest in 
the ellipsoid, and the no-slip condition at the boundary proceeds to set the fluid in rotation. The first peak in W, just 
after t = 0, is due to the Ekman pumping which appears during the spin-up stage, which typically takes place over 
the Ekman time ?£ = E~ l t 2 QT 1 , much faster than the viscous time scale t v = R^/v = E~ l Q _I [3]. For instance, 
in the case shown on figure [2] (a), the dimensionless Ekman time gives Q. t£ ~ 22, which agrees with the numerical 
results. After this initial stage, the fluid is essentially in solid body rotation. From this state, the exponential growth 
of the instability can be seen, before an overshoot and a stationary saturation. Having defined the growth rate <x of the 
elliptical instability as the time constant of the exponential growth, a convergence study on this growth rate is given in 
the figure|2](b). The number of degrees of freedom (DoF) used in most of the simulations of this work ranges between 
4 ■ 10 4 DoF and 7 ■ 10 4 DoF, depending on the ellipticity and the Ekman number, in order to reach a compromise 
between a good convergence (see figure |2](b)) and a reasonable CPU time. 

2.2. Validation of our numerical simulations 

A first visual validation is done on the shape of the flow in comparing the experimental visualisation on figure[3](a) 
with the numerical simulation on figure|3](b). As can be seen, we recover the classical S-shape of the spin-over mode, 
which induces an additional solid body rotation of the fluid around the axis of maximum strain, i.e. perp endicular to 
the imposed rotation axis (Oz) and at an angle of about 45° compared to the deformation axis (Ox) B27I1 . In order to 
quantitatively validate the numerical model, the evolution of the growth rate of the instability is compared in figure |4] 
to the linear theory given in 12711 for small ellipticities : 

°- = \-K<S., (3) 

e 2 e 

where K is a constant equal to K — 2.62 in the limit of small s ( lfl9[|27lo . Note that the second term on the right hand 
side of the equation[3]corresponds to the viscous damping of the growth rate due to the presence of Ekman layers near 
boundaries [see |25l[l9|]. The expression 0) allows to define a critical Ekman number for the onset of the instability 
equal to 



[hi 



(4) 



Close to the threshold, the numerical results closely follow the linear analysis, for values of ellipticity as large as 
s » 0.5. Note in particular that all curves for various (e, E) superimpose providing that cr/e is expressed as a function 
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of yE/s. This is especially interesting in order to apply our results to planets, whose very small E are not directly 
accessible by our numerical tool, but can be compensated for by large s in simulations. 

Finally, for very large ellipticities (e > 0.5), figure |4] shows that the growth rate decreases toward zero. The ellip- 
tical instability finally disappears for ellipticities greater than a critical value which depends on the Ekman number. 
The position of the maximum for the variation of cr with the ellipticity is around s = 0.5, whatever the Ekman number 
is. This indicates that the decrease is probably due to geometrical effects related to the large value of s, rather than to 
any viscous attenuation. 



3. Systematic numerical studies of geophysical complexities 

3.1. Influence of the length of the ellipsoid along the rotation axis 

Having validated our numerical code, we can investigate the influence of the aspect ratio c/b on the instability for 
a given ellipticity in the equatorial plane and a given rotation rate around the axis (Oz). This is directly related to geo- 
and astro-physical flows as oblateness of planets, like the Earth for instance, is most of the time much stronger than 
its tidal deformation, meaning that the rotation axis is also the smallest one. The situation is even more pronounced in 
certain stars, as for instance Regulus A whose diameter is about 32% greater at the equator than the distance between 



its poles [44]. 



Early theoretical work on this aspect was done by [21], who considered the inertial wave basis of an oblate 
spheroid (a = b) and calculated the first 60 subharmonic exact resonances and their growth rate for small ellipticities 
depending on the value of c [see also|23, for a special application to the case of Io]. An explicit theoretical answer for 
the growth rate has also been found by 11511 . starting from the base flow 

a b ,cs 

Ub = -7 ye x + - xt y , (5) 
b a 

and looking for inviscid perturbations that are linear in space variables, corresponding to the classical spin-over mode. 
Here, e x and e y are respectively the unit vectors of {Ox) and (Oy). In an open domain, the base flow © corresponds 
to elliptical streamlines with an ellipticity s = (b 2 - a 2 )/(b 2 + a 2 ) as in our numerical model, but with a variable 
tangential velocity. As shown in figure [5] this base flow is a very good approximation of our configuration where we 
impose a constant tangential velocity on the elliptical boundary, outside a small boundary layer close to the external 
wall where recirculation cells take place. For such a flow, the inviscid growth rate determined by lfl5l1 is: 



Ub 2 - c 2 )(c 2 - a 2 ) 

a \ (b 2 + c 2 )(a 2 + c 2 ) 1 ' 

Note that this theoretical growth rate is valid for a < c < b only and is zero for c = b or c — a. Then the maximal 
theoretical growth rate <x max = ^rf is obtained for c = Vab. One can notice that expression (0 leading to cr = | in 
the inviscid case is recovered with a = R eq + s and b — R eq - s in the limit s — > 0. Note also that the experimental 
choice c — (a + b)/2, equivalent to c — Vab for small ellipticities, results in the growth rate remaining very close to 
the maximum; for instance, even for an ellipticity of 0.7, the difference between the growth rate calculated for c = S±S 
and the maximum value is only 2%. The numerical results are presented in figure[6] In order to compare them with 
the inviscid analytical prediction ©, a viscous damping term —K Ve is added to the expression of the inviscid growth 
rate, similarly to the classical expression (0. Excellent agreement is then found for a < c < b, again validating our 
approach. However, a slightly different constant K = 2.5 (instead of 2.62 in ([§}, valid in the limit of small e) allows a 
better fit of our data. For strong deformations, this constant clearly depends on the considered ellipticity; for instance, 
a similar study performed at s — 0.42 implies K = 2.78. Nevertheless, one can notice that K always remains about 
the same order of magnitude. One can also notice that because of the scattering in the numerical results presented in 
figure|4](a), such a constant K = 2.5 works equally well. Hence in the following, we systematically use K = 2.5. 



In addition to the verification of the theoretical law of 01511 . we are now in the position of exploring the range 
outside a < c < b, where other modes may grow that are not necessarily linear in space variables. As shown in figure 
[6] different modes of the elliptical instability, characterized by their main frequency of oscillation, appear depending 
on the ratio c/b. In this view, the variation of the oblateness is similar to the variation of the aspect ratio in the case of an 
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elliptically deformed cylinder it 1 1 fl . Since the elliptical instability comes from the parametric resonance of two inertial 
waves of azimuthal wave number m and m + 2 with the underlying strain field, the corresponding mode is written 
(to, to + 2) and is related to a pulsation of frequency to mo d e = m + 1 . For instance, the spin-over mode corresponds 
to the stationary mode (—1,1) with half a wavelength along the axis of rotation. According to figure 0(b), the mode 
(1,3) can be observed when c < b and the mode (0,2) when c > a. For even larger aspect ratio, the mode (-1,1) with 
a larger wavenumber than the spin-over takes place. Note however that because of the geometrical confinement, no 
stationary mode can be excited for c < b. 



3.2. Influence of a background rotation 

In geo- and astrophysics, the tidal deformation of a given body (planet, moon or star) is also in rotation because 
of the orbital motion of the companion body. The influence of this background rotation on the development of the 
elliptical instability has been studied theoretically, using short-wavelength analysis 1 8, 33, 35] or normal mode analysis 



Il5il2lll . and also numerically for specific vortices such as Stuart vortices j34ij 5 || or Taylor-Green vortices lf57ll . It 
has been studied experimentally in deformed cylinders i60l [32ll and ellipsoids ikl blll . All these works show that the 
background rotation has a stabilizing effect on cyclones and a destabilizing effect on anticyclones, except when the 
background rotation almost compensates for the flow rotation, in which case the elliptical instability disappears. 

The theoretical expression of the growth rate for perturbations that are linear in space variables (i.e. corresponding 
to the classical spin-over mode) can be readily obtained by following the same method as [15], but taking into account 
an additional Coriolis force coming from the background rotation at a given angular velocity Q c = Q*Q: 



(b 2 - c 2 + 2 Q* a b)(c 2 - a 2 -2Q* a b) 
°~~ \ (b 2 + c 2 )(a 2 + c 2 ) ' (?) 

Actually, this is a particular case of the more general stability analysis detailled in f3] for a precessing triaxial ellipsoid. 
The existence range of this mode is now between c - b Jl + 2D.*a/b and c - a Jl + 2 Q* b/a and the maximal 
theoretical growth rate is now exactly obtained for ^-r — Q* + -^/l + Q* (a/b + b/a + Q*). In our numerical model, 
the addition of such a Coriolis force is straightforward and numerical results are presented in figure[7](a). Again, they 
compare well with the theoretical prediction, providing that the viscous damping term, -2.5 Ve, found in section [37X1 
is added to the inviscid expression (|7). Other modes are selected by the Coriolis force outside of the spin-over mode 
resonance band, as already observed experimentally by ll3lll . 

Finally, we can investigate how the oblateness and the background rotation interact with each others. Considering 
for instance an oblate ellipsoid with a ratio c/b such as the mode (1, 3) is selected in absence of background rotation, 
the spin-over mode is recovered when decreasing Q c , as shown in figure[7](b). More generally, we find that the (—1,1) 
mode is the most unstable one in the anticyclonic domain up to a value Q. c /Q. w , < — 1 where the flow restabilises, in 



agreement with the conclusions of [ 3 1 ] 



3.3. Influence of the obliquity of the ellipsoid 

In planetary cores, tidal deformations are aligned with the orbital plane rather than with the equatorial plane, which 
means that they are not orthogonal to the rotation axis. Previous works do not take into account this phenomenon, 
even though the obliquity can be significant, e.g. 23°26' for the Earth. In this section, we thus investigate the effect of 
obliquity for the first time and consider that, as a first approximation, the shape of the body remains a triaxial ellipsoid. 

In the numerical model, the rotation axis is tilted and oriented along the unit vector 

k c = (cos(0) sin(6i), sin(0) sin(0), cos(6>)), 

where <p is its azimuth angle and 9 its colatitude angle in spherical coordinates (fig. [TJ. In assuming for example that 
the rotation axis is tilted in the (Oxz) plane (i.e. <p = 0), we can study how the development of the elliptical instability 
changes depending on the obliquity 9. For 9 = 0, one recovers the usual configuration with s = (b 2 -a 2 ) / (b 2 +a 2 ) and an 
aspect ratio c/b. For 9 = n/2, one recovers results from section !3. 1 l in exchanging a and c. i.e. e = \b 2 -c 2 \/(b 2 +c 2 ) and 
the aspect ratio is now a/b. In between, streamlines in planes perpendicular to the rotation axis are also elliptical, but 
their centers are not located along the rotation axis anymore, except in the equatorial plane z = 0. Besides the effective 
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ellipticity measured in each plane now depends on z, the maximum being reached in the equatorial plane z — 0. Seen 

■ 2 ,2 - 1 

from this equatorial plane, the apparent length of the ellipsoid along the rotation axis is c - (^M + 2 , whereas 

the great and small axes of the elliptical streamlines are respectively b = b and 5 = (=2^-2 + " . One can then 

estimate the growth rate of the instability using formulas (O or (|6]) with these apparent lengths. Numerical results for 
the spin-over mode compared with these approximations are presented in figure [8] Results agree well for obliquity 
up to 6 ~ 20°. Then, the numerical growth rate significantly differs, probably because the geometry is very far from 
an apparent ellipsoid in rotation around one of its principal axes. It remains however between the two expressions 
proposed. Note that for the Ekman number studied here, we do not see the tidal instability reappearing around 9 ~ 90°, 
which will however be the case at smaller Ekman number: as seen in section lXTl in this case other modes could appear 
because of the modified oblateness. 



4. Scaling laws for global quantities of geophysical interest 

4.1. Scaling law for the amplitude of the flow driven by the instability 

The amplitude of the flow driven by the instability at saturation remains up to now an open question, mainly be- 
cause it is determined by strong non-linear interactions. It is however an important quantity for geo- and astrophysical 
applications, since it allows us to evaluate the influence of elliptical aspects compared to the other relevant ingredients 
of planetary dynamics, such as convection. In order to study the amplitude of the flow driven by the instability at 
saturation, we define the amplitude A* by the maximum value over the volume V: 

A* = max ||u - Ubll (8) 

v 



where u is the dimensionless velocity field and Ub is the theoretical base flow defined in section [3TI The evolution 
of A* as a function of the Ekman number is shown in figure [9] (a) for various ellipticities: it can be seen that A* is 
not zero below the threshold of the elliptical instability because of the differences around the outer bound between 
the theoretical base flow and the flow numerically obtained (see figure |5J- Actually, the difference is maximal on 
the boundary, along the smallest equatorial axis where the theoretical velocity is maximal. There, the corresponding 
amplitude A( S ) can be calculated theoretically : 

b 2 
~ 1 + 

which only depends on e. One then defines the amplitude of the flow driven by the instability by A = A* - A( B) , which 
is equal to below the threshold (fig. |9](b)). 

Far from threshold, a secondary instability appears, which induces a secondary dynamics superimposed on the 
primary state : e.g. typically, for s = 0.317, E = 1/1500, c - (a + b)/2, the spinover mode is no more stationary, and 
the flow is slightly oscillating around the spinover mean flow, at a pulsation (±> sec » 1 .4. Note that this is in agreement 
with I22I1 which predicts that the primary elliptical instability should saturate and be stable only in a small strain 
window e - e c — O(E). Thus, the amplitude A has to be averaged in time to smooth out the small scale fluctuations. 
Since those small-scale fluctuations take place on a typical time scale comparable to one revolution In/Q., whereas the 
characteristic time for the tidal instability is the inverse of the growth rate, i.e. (Q e/2) _1 , the averaging is performed 
about a typical time Q. 1 yfAnjE corresponding to the geometrical mean of these two extreme values. 

With this definition, figure [9] (b) shows that all results collapse on the same generic law above the threshold 
providing that we use the variable E c /E - 1. Far from threshold, the amplitude seems to saturate around 1, which 
means that the velocities generated by the tidal instability are comparable to the imposed boundary rotation. 

One can notice that near the threshold, a square root A « 0.6 yjE c IE - 1 fits the results, which is in agreement 
with a pitchfork bifurcation. Actually, this scaling law can be obtained analytically, starting from the simple model 



used in B27H to describe the viscous non-linear evolution of the spinover mode. This model, first introduced by [2 1 



and 114911 for an inviscid solid-body rotation in a spheroid, reduces to: 

co x — —a (1 + co z ) a)y — v so a> x (10) 
6 



(jjy — -/3 (1 + L0 Z ) co x - v so to- 



(11) 



u) z = s o) x (jjy - v ec a) z + v„i (cj 2 x + L0 2 y ) (12) 

where co is the rotation vector of the spinover mode, a = and /3 = ^- . The damping terms are given by theory 
(no adjustable parameter): v so = 2.62 y/E is the linear viscous damping rate of the spinover mode (first calculated 
by il2l0 . v ec = 2.85 y/E is the linear viscous damping of axial rotation and v„i — 1.42 y/~E is the viscous boundary 
layer effect on the non-linear interaction of the spinover mode with itself (see IU2I0 . Even if this model does not 
take into account all the viscous terms of 0( VE) or the non-linear corrections induced by the internal shear layers, 



it satisfyingly agrees with experiments, regarding the growth rate as well as the non-linear saturation [27]. Thus, this 
model can be used with confidence to describe the viscous non-linear evolution of the spinover mode. 
After little algebra, we obtain a non-trivial stationary state for e > 2v so , given by: 



_ + i Vec [ ~ Vso] _ / v ec [s-2 y~J 

Ux ~ ~ V/3 s - v„, [ V^8 +/3 2 / Vo/3] ~ " V s 1 - 2 v nl s 



v ec [y/a]3-v so ] _ v ec [s-2 v so ] _ 

0)y = + -t — " ■=■ ~ + ,/ ~ + Wi (14) 

\ a s- v„i [ yfaf} + a 2 1 y[a]S\ V e - 2 v„i s 

cj z = ^- y/T^-l^^-l (15) 
s s 

where approximations are done assuming s «: 1 . Now, the amplitude A corresponds to the norm of the flow a> x r 
driven by the spinover mode. Then, near the threshold (e m 2 ■ 2.62 V^e and E c /E - 1 « 1), we obtain: 



(16) 



which is in good agreement with the numerical fit. 



4.2. Scaling law for the viscous dissipation by the instability 

As explained for instance in lf5lll or in jjlll . the energy dissipated by tides is a primordial quantity which directly 
influences the orbital evolution and rotational history of a binary system during its synchronization. It is however 
poorly known [e.g. |63l 15911 . In most traditional models, fluid dissipation in the planetary core is supposed to be 
negligible. However this may not be the case when the elliptical instability is excited at a rather large amplitude, 
inducing important shears between the bulk of the fluid and the boundary. Our purpose here is to systematically 
quantify the variation of this dissipation with the ellipticity and the Ekman number. 

The dissipated power balance of the incompressible flow in the ellipsoid is given by: 



at 



(pyjdr= JjT (5yn). ft dj- JJJ «5" v : Vu dr (17) 



where cr v = 77 (Vu + ' Vu) is the viscous stress tensor of the newtonian fluid, S the surface of the ellipsoid, u = Q R eq u 
is the dimensionalized velocity field, p the volumic mass and 77 the dynamic viscosity of the fluid. In this section, we 
focus on the stationary spinover mode and then the two terms on the right side of equation ( TTTb balance each others: 
the first one allows to maintain the rotation of the fluid by the no-slip conditions on the boundary while the second 
one is the volume dissipation P of the fluid. In the following, we use as a power scale the dissipated power during the 
spin-down stage of the equivalent sphere of radius R eq and moment of inertia 1^(1^-2/5 MR 2 for an homogeneous 
sphere of mass M and radius R), i.e. the kinetic energy of rotation E c — 1 Q 2 divided by the Ekman time f£. Note 



that most of the dissipated power comes naturally from the boundary layers: its determination is thus a challenging 
task from a numerical point of view, and extra care must be taken regarding the convergence of the simulations, as 
shown for instance in figure [2] (b). 

Below the threshold, the dissipated power is not zero because of the recirculation patterns of our base flow (see 
figure[3). This dissipation is due to the flow driven by the ellipticity, which scales as s OR for small ellipticities. Then, 
the dimensionalized dissipated power below the threshold simply scales as 

<x v : Vu dr ~ rj s 2 CI 2 R 3 (18) 



v 



which is confirmed by the dissipated power measured in our numerical simulations below the threshold. 



Far from threshold, the model proposed by [31] considers that the spin-over mode simply corresponds to a sup- 



plementary solid body rotation Q§o outside the outer viscous boundary layer of thickness h — £ vv/Aso, £ being a 
constant of order 1 . Then, the dissipation is only located in this boundary layer, where the fluid rotation has to match 
the imposed velocity conditions at the outer boundary. With this simple model, the torque of the container on the fluid 
is C m / C = — 2 M v f ^so> an d the power dissipated by the system is P — -2 M v | £lso 2 , which can also be written 
P — — 1 M R sfv £l 5/2 A 5 ^ 2 , where A is the dimensionless spin-over mode amplitude studied in the previous section. 
The dimensionless power given by this model can finally be stated as 

p = \£\^ 1 1 A V2 (19) 

P \ I A Q 3 V£ f 

where, according to the previous section, A is about 1 in the small Ekman number limit studied here. 

The evolution of Pdissip measured in our numerical model is shown in figure [10] It confirms the behavior pre- 
dicted by (TT~9b with a saturation far from threshold and ^ ~ 1 . This behavior is especially interesting for planetary 
applications, as will be studied in the following section. 



5. Orders of magnitude for planetary applications 

From a geophysical point of view, the results from the previous sections can be used to derive the orders of 
magnitude involved in planetary cores. Regarding the jovian moon Io, which is clearly unstable to the tidal instability 
and where the strong tidal dissipation is a topical question 1 30], our numerical study confirms the first trends given by 
lUlll . with however a numerically determined scaling factor of the order l/£ a 1 in the power dissipation estimates. 
Moreover, the most important result comes from the oblateness effects. Indeed, according to the data given in ll23ll . 
c/b ~ 0.995 and a/b > 0.999, which means that the excited mode of the tidal instability in Io cannot be the spin- 
over mode, but has to be oscillatory. Note that this result is not modified when the background rotation is taken into 
account. 

A similar analysis can be done for the early Earth, with a Moon two times closer than today. In this case, the 
length of the day was around 10 hours according to lf59ll and the lunar tides were around eight times stronger, whereas 
the solar tides were about the same. Then, the actual tidal amplitude of 50 cm allows to estimate the ellipticity of 
the early Earth: s ~ IO -6 . In considering a similar but totally molten core, we estimate E c /E - 1 ~ 70, which 
means that the early Earth's core was clearly unstable to tidal instability, and that the tidal instability was then at 
saturation. Moreover, with an orbital period for the Moon (1 /2) 3 ^ 2 shorter than the actual moon (Kepler's Third Law), 
the dissipated power given by the model was around P ~ 10 18 W. This estimation seems huge in comparison with 

1 7 J li 

the present dissipation by tidal friction (~ 3.75 • 10 W according to [46]) but actually, it simply suggests that the 
Earth-Moon system was then in rapid evolution. Once again, the possible mode of the tidal instability was not the 
spin-over mode neither any stationary (—1,1) mode, considering that the oblateness of the rapidly rotating early Earth 
was larger than the actual oblateness. In fact, we expect stationary modes to be only marginally excited in geophysical 
systems, limited to very peculiar configuration regarding the oblateness and the rotation of the tidal bulge. 

Note that computations similar to the ones presented here for an ellipsoidal geometry have been performed for 
an ellipsoidal shell in order to study the influence of a solid inner planetary core. Similar conclusions are then found 
regarding the physics of the tidal instability. Only a slightly larger dissipation term due to the presence of an inner 
viscous boundary layer has to be taken into account. 
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6. Conclusion 



This paper presents the first systematic numerical study of the tidal instability in a rotating ellipsoid. The numerical 
approach is a powerful tool to complete the knowledge derived from previous theoretical and experimental studies. 
The effects of oblateness, background rotation and obliquity, which lead to the selection of various instability modes, 
have been studied. We have also defined scaling laws regarding the amplitude of the flow driven by the instability and 
its viscous dissipation, which confirm the primordial role played by tidal effects at a planetary scale. 

Another fundamental geophysical issue regarding the elliptical instability is whether or not it can drive a planetary 
dynamo. The results presented herein show that the numerical approach is able to faithfully capture the physics of 
the hydrodynamic elliptical instability. Thus, the present work provides a solid basis to study the full MHD problem 
encountered in planetary core flows undergoing tidal instability. These flows can only be studied numerically due to 
their inherent complexity. Therefore, a next step will be to introduce the effects of a magnetic field in our model. 
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Figure 2: Simulations performed at E = 1/500 and e = 0.317 for c = (a) Time evolution of the mean value of the vertical velocity, showing the 
spin-up phase (CI tg as 22), the exponential growth of the tidal instability (see also in dotted line the exponential fit with a growth rate cr = 0.352) 
and its saturation, (b) Convergence with the number of degrees of freedom (DoF) of the growth rate cr and of the dissipated power Pdissip at 
saturation. Results presented in (a) are computed for 42459 DoF. 
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(b) 

Figure 3: Validation of the numerical simulations, (a) Kalliroscopic visualization of the spin-over mode in the meridional plane of maximum shear 
for E = 1/4000 and e = 0.16. The typical S shape of the rotation axis is due to the combination of the main rotation imposed by the boundary 
and the spin-over mode, (b) Slices of the velocity field ||u|| and surface of iso-value ||u|| = 0.15 at saturation of the tidal instability for E = 1/344, 
e = 0.317, c = St^j and 49900 DoF. The classical S-shape of the spin-over mode is recovered. 
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Figure 4: Validation of the numerical simulations, (a) Evolution of the growth rate for (E > 1/2000, £ < 0.95, c = ^ir) and comparison with 
the linear theory indicated by a dashed line. The coefficient 2.62 comes from KFhl and is valid in the limit of small ellipticity. Good agreement 
is found for values of e up to 0.5. (h) Evolution of the growth rate depending on the ellipticity for two values of the Ekman number E = 1/344 
and E = 1/800 (e < 0.8, c = As also seen in (a), the growth rate agrees with the linear analytical analysis close to the threshold and then 
decreases for large values of the ellipticity. 
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Figure 5: Azimuthal velocity vg for E = 1/95, e = 0.1, c = 2±£. (a) Slice along (Oy) at x = 0. (b) Slice along (Ox) at y = 0. (c) Slice in 
the equatorial plane of the vertical component of the vorticity. Dashed lines in (a) and (b) correspond to the theoretical base flow 0, which 
presents a variable tangential velocity along elliptical streamlines. Good agreement is found with the numerical results, except in the small outer 
viscous boundary layer, where recirculation cells take place to match the imposed constant velocity along the boundary. Note in particular that the 
maximum velocity is not reached at the boundary but within this viscous boundary layer. 
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Figure 6: Influence of the length of the ellipsoid along the rotation axis on the mode selection of the tidal instability (E = 1/688, £ = 0.317J. (a) 
Variation of the growth rate, (h) Variation of the main frequency of the selected mode determined by Fourier analysis of its saturation state. 
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Figure 7: Evolution of the growth rate in the presence of a background rotation for a fixed ellipticity e = 0.317 and a fixed value of the total Ekman 
number E to , = — y -^r- = 10~ 3 , where Sl m takes into account the background rotation and the fluid rotation, i.e. fl fof = £l c + fl. (a) c = (a + b)/2, 

O-tot %eq 

i.e. c/b = 0.86: the spin-over mode is then excited in the absence of background rotation. Good agreement is found around this value with the 
analytical solution (7). Further decreasing Sl c , other modes with smaller wavelength along the rotation axis appear, (b) c/b = 0.65: the (1,3) 
mode is then excited in the absence of background rotation. Other modes can be excited: • represent the (1,3) mode, A the spin-over mode, ■ 
represent the (—1,1) mode with one wavelength along the axis of rotation, and * are other modes. 
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(a) 

Figure 8: Evolution of the growth rate with the obliquity 8 for E = 1/600 (fixed angular rate) and two values of the ellipticity e = 0.317 and 
£ = 0.42. The dashed dotted line and the dashed line correspond to expression J§), taking into account the values of apparent axes seen from 
the equatorial plane, whereas the dotted line and the continuous line correspond to expression 0, simply considering the apparent value of the 
ellipticity seen from the equatorial plane. The coefficient Kfor viscous corrections is determined at0 = O° and then kept constant. 
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Figure 9: Variation of the saturation amplitude of the flow driven by the elliptical instability depending on the Ekman number, for various values 
of the ellipticity. (a) Maximum value of the difference between the actual velocity and the theoretical base flow 0. (b) Maximum amplitude of 
the tidal instability, corresponding to the previous values corrected to take into account recirculation cells: all curves then superimpose when 
computed as a function of'E c /E- 1. 
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Figure 10: Viscous dissipation by the tidal instability, as a function of the distance from threshold. The results collapse on a generic law and seem 
to converge towards a saturation value far from threshold. 
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